###
# Julia script for replicating quantitative model analysis in "Q-Monetary Transmission", by Priit Jeenas (priit.jeenas@upf.edu) and Ricardo Lagos (ricardo.lagos@nyu.edu). May 30, 2023.
###
# Employ homogeneity of degree 1 in incoming capital, and solve over outstanding equity (s) space: Ke (s_{t+1}=S_{t+1}/K_{t+1}) as choice for next period
# NOTE: dir_path must be the local directory of the replication package
dir_path = "SET_TO_LOCAL_DIRECTORY"

cd(dir_path)
push!(LOAD_PATH, pwd())

# Load procedures written for solving the model
using QMT_ihzn_1_xirs
# Load packages
using Plots
Plots.gr()
using LaTeXStrings
using Roots
using CSV
using Distributions
using DataFrames

# Setting up parameters
# Calibration (exogenously set parameters)
beta_set    = 0.995
eta_set     = 2.0   # Exponent on convex capital adjustment cost (2.0 = quadratic)
delta_set   = 0.025 # Depreciation rate
delta_0_set = 0.025 # "Center point" of convex capital adjustment cost
pi_set      = 1.0-0.017 # Quarterly survival rate
rm_set      = (1.0/beta_set-1.0) # Real return on firms' "liquid assets" (m) held as bonds
m_0_set     = 0.40/0.60 # Entrants' (IPO)

# Calibrated parameter values
eps_e_set, nu_set, rK_set, xi_ub_set, a_set = 4.0083, 45.3184, 0.0195, 0.145,  0.0289 # Entrepreneur's valuation of dividends (eps_e), scalar multiplying convex adjustment costs (nu), capital productivity in subperiod 1 (rK), upper bound on equity issuance costs uniform distribution (xi_ub), capital productivity in subperiod 2 in the quantitative extension (a)

# Code allows to introduce "shock" to firms' liquid asset holdings - switch off probability to zero
prob_msh_set= 0.0
x_msh_set   = 0.0

# Determine psi_set (steady state stock price) based on the Euler equations
if true
    iota_ss = 0.04/4            # Steady state quarterly nominal rate
    sigma_eps   = 2.56          # sigma in the outside investors' LogNormal dividend valuation shock distribution
    mu_eps      = -0.5*sigma_eps^2 # mu in the outside investors' LogNormal dividend valuation shock distribution
    eps_bar = exp(mu_eps+0.5*sigma_eps^2)
    alpha   = 1.0               # Normalization
    theta   = 1.0               # Normalization
    # Initiate standard normal
    snd     = Normal()

    # SS value of psi, as function of eps_star
    function fun_psi_ss(eps_star)
        return (beta_set*rK_set/(1-beta_set*pi_set*(1-delta_set)))  * (pi_set*a_set/rK_set + eps_bar + alpha*theta*(eps_star*cdf(snd, (log(eps_star)-mu_eps)/sigma_eps) - eps_bar*cdf(snd, (log(eps_star)-(mu_eps+sigma_eps^2))/sigma_eps) ) )
    end

    # Analyze RHS of SS EEm as function of eps_star (x) and psi_next (y). This is a general function that holds outside of SS also
    rhs_EEm(x,y) = alpha*theta*(1.0/(x + pi_set*((1-delta_set)*y + a_set)/rK_set)) * (eps_bar*cdf(snd, (mu_eps+sigma_eps^2-log(x))/sigma_eps) -x*(1.0-cdf(snd, (log(x)-mu_eps)/sigma_eps)))
    # Because in the range (eps_L, eps_H), the RHS is strictly decreasing in eps_star, there is a unique solution to ε*∈[ε_L,ε_H] as long as ι is between the following bounds
    iota_L  = rhs_EEm(1e8, fun_psi_ss(1e8))
    iota_H  = rhs_EEm(0.0, fun_psi_ss(0.0))

    # Function to solve for eps_star as a function of iota_ss
    function fun_eps_star(iota_val)
        return find_zero(x -> iota_val - rhs_EEm(x, fun_psi_ss(x)), (0.0, 1e8), Bisection() )
    end
    # Method that takes next period psi as given
    function fun_eps_star(iota_val, psi_n)
        return find_zero(x -> iota_val - rhs_EEm(x, psi_n), (0.0, 1e8), Bisection() )
    end

    Plots.plot(x -> rhs_EEm(x, fun_psi_ss(x)), 0.0, 10.0)
    Plots.plot(x -> fun_eps_star(x), 1e-7, iota_H, label="")

    # Solve for the actual values
    eps_star_ss = fun_eps_star(iota_ss)
    psi_ss      = fun_psi_ss(eps_star_ss)
    # Impose in model
    psi_set     = psi_ss
end

# Set bounds on grids:
m_lb, m_ub  = 0.0, 0.7      # Bounds on liquidity ratio grid
i_ub        = 0.2           # Upper bound on investment rate choice grid
Kep_ub      = 1-delta_set+i_ub   # Natural upper bound on Kep choice (=s_{t+1}/k_{t+1}) (equity outstanding relative to k going forward cannot exceed 1)

if true
    spliorder_set = [1,1]
    n_set       = [80,5]
    nf_set      = [10,100,40]
    np_set      = [80,20,80]
    curv_set    = [0.5,1.0]
    curvf_set   = [1.0,1.0,1.0]
    curvp_set   = [1.0,1.0,1.0]
    curvp_set   = [0.7,0.7,0.7]

    mmin_set, mmax_set  = m_lb, m_ub
    imin_set, imax_set  = 0.0, i_ub
    ymin_set, ymax_set  = 0.0, m_ub
    Kepmin_set, Kepmax_set = 0.0, Kep_ub
    Kemin_set, Kemax_set   = 0.0, 1.0
    Kminf_set, Kmaxf_set= 0.1, 5.0
    mminf_set, mmaxf_set= m_lb, m_ub
    Keminf_set, Kemaxf_set= 0.0, 1.0
end

# Create empty text file to save model output in "mod_savedout" folder
if !isdir("mod_savedout") mkdir("mod_savedout") end
if !isdir("mod_savedout/conv_figs") mkdir("mod_savedout/conv_figs") end
modlname = "mod"
ofile = open("mod_savedout/"*modlname*"_output.txt", "w")
close(ofile)

# Construct solution parameter container "opts" and state space container "sspace"
opts = Options(Nbell=200, Nhow=100, first_how=10, tolc=1e-7, itermaxL=1000, tolL=1e-8, plotSD="N", plotc="N", sfigs="Y");
sspace = StateSpace(spliorder=spliorder_set, n=n_set, nf=nf_set, np=np_set, curv=curv_set, curvf=curvf_set, curvp=curvp_set, mmin=mmin_set, mmax=mmax_set, Kemin=Kemin_set, Kemax=Kemax_set, imin=imin_set, imax=imax_set, ymin=ymin_set, ymax=ymax_set, Kepmin=Kepmin_set, Kepmax=Kepmax_set, Kminf=Kminf_set, Kmaxf=Kmaxf_set, mminf=mminf_set, mmaxf=mmaxf_set, Keminf=Keminf_set, Kemaxf=Kemaxf_set);
modl = Model(beta=beta_set, eps_e=eps_e_set, rK=rK_set, a=a_set, nu=nu_set, eta=eta_set, delta=delta_set, delta_0=delta_0_set, pi=pi_set, psi=psi_set, rm=rm_set, m_0=m_0_set, prob_msh=prob_msh_set, x_msh=x_msh_set, xi_ub=xi_ub_set, name=modlname);

# Setup model parameters and state space containers
modl, sspace = setup!(modl, sspace, opts);

# Solve firm problem in steady state with loop
c_sol = solve_c(modl, sspace, opts, [], [])

# Compute steady state distribution of firms across the liquidity ratio and the equity outstanding (relative to K) ratios
opts.plotSD="Y"
L_sol = solve_L(c_sol, modl, sspace, opts)

####
# Compute IRF to MIT shock to nominal rate (iota) path
####
using Distributions
using Statistics

# Set seed for random number generator in sampling firms
import Random
Random.seed!(1234)
# Number of firms and periods
Nf          = 20000
T           = 13
manl_ftsize = 13

# Compute SS policies
v_sol1, = solve_valfunc_v(c_sol, modl.psi, modl.rm, sspace.s, modl, sspace, opts, false)
v_sol0, = solve_valfunc_v0(c_sol, modl.psi, modl.rm, sspace.s, modl, sspace, opts, false)
ss_xi_vec = v_sol1.val-v_sol0.val
# And interpolation coefficients for policies in steady state
i_cpol1, mp_cpol1, Kep_cpol1, ie_cpol1, y_cpol1 = sspace.Phi_v \ v_sol1.i, sspace.Phi_v \ v_sol1.mp, sspace.Phi_v \ v_sol1.Kep, sspace.Phi_v \ v_sol1.ie, sspace.Phi_v \ v_sol1.y
i_cpol0, mp_cpol0, Kep_cpol0, ie_cpol0, y_cpol0 = sspace.Phi_v \ v_sol0.i, sspace.Phi_v \ v_sol0.mp, sspace.Phi_v \ v_sol0.Kep, sspace.Phi_v \ v_sol0.ie, sspace.Phi_v \ v_sol0.y
ss_xi_cvec = sspace.Phi_v \ ss_xi_vec

# Collections of shocks to average across (Write them as two shocks of +/- 25bp, but below adjust paths to imply a +/- 1% stock price response)
bp_shocks_vec = [25.0, -25.0]
sweights_vec  = [0.50, 0.50] # Weights to give shocks in final averaging across scenarios
sweights_vec  = sweights_vec / sum(sweights_vec)

Nshocks  = length(bp_shocks_vec)
# Create collectors of MEAN IRFs, for [all, high, low]. For each shock value.
firms_i_mean_IRFdev_col, firms_mp_mean_IRFdev_col, firms_ie_mean_IRFdev_col, firms_y_mean_IRFdev_col, firms_Ei_mean_IRFdev_col, firms_li_mean_IRFdev_col  = zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T)
firms_i_mean_IRFdev_high_col, firms_mp_mean_IRFdev_high_col, firms_ie_mean_IRFdev_high_col, firms_y_mean_IRFdev_high_col, firms_Ei_mean_IRFdev_high_col, firms_li_mean_IRFdev_high_col  = zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T)
firms_i_mean_IRFdev_low_col, firms_mp_mean_IRFdev_low_col, firms_ie_mean_IRFdev_low_col, firms_y_mean_IRFdev_low_col, firms_Ei_mean_IRFdev_low_col, firms_li_mean_IRFdev_low_col  = zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T), zeros(Nshocks, T)
# Also for log AC
firms_liAC_mean_IRFdev_col  = zeros(Nshocks, T)
firms_liAC_mean_IRFdev_high_col  = zeros(Nshocks, T)
firms_liAC_mean_IRFdev_low_col   = zeros(Nshocks, T)

### Start loop
for is=1:Nshocks
    bp_shock    = bp_shocks_vec[is]

    println("Running IRF for bp_shock=$(bp_shock).")

    # Compute IRF
    if true
        # Compute psi_col based on nominal rate shock
        rho_iota = 0.5  # Persistence of psi_shock

        # Scaling manually for 25bp shocks, so that expansionary shock leads to 1% increase in stock price, contractionary to 1% decrease
        if is==1
            eps_scale = 1.19 # Use this "eps_scale" here to "pre-scale" eps_iota for getting a 1% psi response for model name 0729-01 (29 July 2022)
        elseif is==2
            eps_scale = 1.0 # Use this "eps_scale" here to "pre-scale" eps_iota for getting a 1% psi response for model name 0729-01 (29 July 2022)
        end
        eps_iota = -(1.0/eps_scale)*bp_shock/(4*10000) # bp change in annualized iota when rho=0.75 (to get 1% stock price change)
        # Collection of iota_{t+1}
        iota_col = iota_ss*ones(T) + [eps_iota*rho_iota^(tt-1) for tt=1:T]
        # Iterate backwards to solve for psi_col
        function fungen_psi_col()
            psi_col  = zeros(T)
            psi_next = psi_ss
            for tt=T:-1:1
                eps_star_next = fun_eps_star(iota_col[tt], psi_next)
                psi_cur       = modl.beta*(eps_bar*modl.rK + modl.pi*modl.a + modl.rK*alpha*theta*(eps_star_next*cdf(snd, (log(eps_star_next)-mu_eps)/sigma_eps) - eps_bar*cdf(snd, (log(eps_star_next)-(mu_eps+sigma_eps^2))/sigma_eps) ) + modl.pi*(1.0-modl.delta)*psi_next) # -- LOGNORMAL ε
                psi_col[tt] = psi_cur
                psi_next    = psi_cur
            end
            return psi_col
        end
        psi_col = fungen_psi_col()
        # Percentage deviations in stock prices
        psi_pct_dev = 100*psi_col./(psi_ss*ones(T)) - 100*ones(T)
        eps_scale = psi_pct_dev[1] # Adjust slightly the final scaling of impulse responses to correspond more precisely to a 1% stock price change
        manprintln("eps_scale=$(eps_scale)", "mod_savedout/"*modl.name*"_output.txt")
    end

    # Compute IRF solution
    irfsoln= solve_IRF(psi_col, c_sol, modl, sspace, opts, [], [])

    # Create interpolation coefficients for policies
    # IRF
    i_cpol1_irf, mp_cpol1_irf, Kep_cpol1_irf, ie_cpol1_irf, y_cpol1_irf = sspace.Phi_v \ irfsoln.i_col1, sspace.Phi_v \ irfsoln.mp_col1, irfsoln.Kep_col1, irfsoln.ie_col1, irfsoln.y_col1
    i_cpol0_irf, mp_cpol0_irf, Kep_cpol0_irf, ie_cpol0_irf, y_cpol0_irf = sspace.Phi_v \ irfsoln.i_col0, sspace.Phi_v \ irfsoln.mp_col0, irfsoln.Kep_col0, irfsoln.ie_col0, irfsoln.y_col0
    irf_xi_cvec_col = sspace.Phi_v \ irfsoln.xi_col

    # Repeat for SS distribution, and split based on median
    si_dist_ss  = wsample(1:sspace.Nsf0, L_sol, Nf)
    s_dist_ss   = sspace.sf0[si_dist_ss,:]
    m_cut_ss    = quantile(s_dist_ss[:,1], 0.5)
    # Draw issuance costs also
    xi_draws_samp = rand(Uniform(0,modl.xi_ub), (Nf,size(irfsoln.i_col1,2)))

    s_dist = s_dist_ss
    Nsamp = size(s_dist,1)

    # Collectors
    firms_i_col_SS, firms_mp_col_SS, firms_ie_col_SS, firms_y_col_SS, firms_Ei_col_SS = zeros(Nsamp, T),  zeros(Nsamp, T), zeros(Nsamp, T), zeros(Nsamp, T),  zeros(Nsamp, T)
    firms_i_col_IRF, firms_mp_col_IRF, firms_ie_col_IRF, firms_y_col_IRF, firms_Ei_col_IRF = zeros(Nsamp, T),  zeros(Nsamp, T), zeros(Nsamp, T), zeros(Nsamp, T),  zeros(Nsamp, T)
    firms_iAC_col_SS = zeros(Nsamp, T)
    firms_iAC_col_IRF = zeros(Nsamp, T)

    # Simulate each single firm
    for ii=1:Nsamp
        s_init_temp     = s_dist[ii,:]'
        xi_draws_temp   = xi_draws_samp[ii,:]
        firm_temp_sim_SS = sim_firm_pols_intp(s_init_temp, T, xi_draws_temp, repeat(i_cpol1, 1, T), repeat(mp_cpol1, 1, T), repeat(Kep_cpol1, 1, T), repeat(ie_cpol1, 1, T), repeat(y_cpol1, 1, T), repeat(i_cpol0, 1, T), repeat(mp_cpol0, 1, T), repeat(Kep_cpol0, 1, T), repeat(ie_cpol0, 1, T), repeat(y_cpol0, 1, T), repeat(ss_xi_cvec, 1, T), modl, sspace, opts)
        firm_temp_sim_IRF = sim_firm_pols_intp(s_init_temp, T, xi_draws_temp, i_cpol1_irf, mp_cpol1_irf, Kep_cpol1_irf, ie_cpol1_irf, y_cpol1_irf, i_cpol0_irf, mp_cpol0_irf, Kep_cpol0_irf, ie_cpol0_irf, y_cpol0_irf, irf_xi_cvec_col, modl, sspace, opts)

        firms_i_col_SS[ii,:], firms_mp_col_SS[ii,:], firms_ie_col_SS[ii,:], firms_y_col_SS[ii,:], firms_Ei_col_SS[ii,:] = firm_temp_sim_SS.i_col, firm_temp_sim_SS.mp_col, firm_temp_sim_SS.ie_col , firm_temp_sim_SS.y_col, firm_temp_sim_SS.Ei_col
        firms_i_col_IRF[ii,:], firms_mp_col_IRF[ii,:], firms_ie_col_IRF[ii,:], firms_y_col_IRF[ii,:], firms_Ei_col_IRF[ii,:] = firm_temp_sim_IRF.i_col, firm_temp_sim_IRF.mp_col, firm_temp_sim_IRF.ie_col , firm_temp_sim_IRF.y_col, firm_temp_sim_IRF.Ei_col

        firms_iAC_col_SS[ii,:]    = firms_i_col_SS[ii,:] + modl.nu * (firms_i_col_SS[ii,:]-modl.delta_0*ones(size(firms_i_col_SS[ii,:],1))).^modl.eta
        firms_iAC_col_IRF[ii,:]   = firms_i_col_IRF[ii,:] + modl.nu * (firms_i_col_IRF[ii,:]-modl.delta_0*ones(size(firms_i_col_IRF[ii,:],1))).^modl.eta
    end

    # Compute means across firms
    lows_vec  = s_dist_ss[:,1] .<= m_cut_ss
    highs_vec = s_dist_ss[:,1] .> m_cut_ss
    firms_i_mean_SS, firms_mp_mean_SS, firms_ie_mean_SS, firms_y_mean_SS, firms_Ei_mean_SS  = vec(mean(firms_i_col_SS, dims=1)), vec(mean(firms_mp_col_SS, dims=1)), vec(mean(firms_ie_col_SS, dims=1)), vec(mean(firms_y_col_SS, dims=1)), vec(mean(firms_Ei_col_SS, dims=1))
    firms_i_mean_IRF, firms_mp_mean_IRF, firms_ie_mean_IRF, firms_y_mean_IRF, firms_Ei_mean_IRF = vec(mean(firms_i_col_IRF, dims=1)), vec(mean(firms_mp_col_IRF, dims=1)), vec(mean(firms_ie_col_IRF, dims=1)), vec(mean(firms_y_col_IRF, dims=1)), vec(mean(firms_Ei_col_IRF, dims=1))
    # Also for log investment rates
    firms_li_mean_SS = vec(mean(log.(firms_i_col_SS), dims=1))
    firms_li_mean_IRF = vec(mean(log.(firms_i_col_IRF), dims=1))
    # Also for AC
    firms_liAC_mean_SS  = vec(mean(log.(firms_iAC_col_SS), dims=1))
    firms_liAC_mean_IRF = vec(mean(log.(firms_iAC_col_IRF), dims=1))

    # Repeat for low subsample
    firms_i_mean_low_SS, firms_mp_mean_low_SS, firms_ie_mean_low_SS, firms_y_mean_low_SS, firms_Ei_mean_low_SS  = vec(mean(firms_i_col_SS[lows_vec,:], dims=1)), vec(mean(firms_mp_col_SS[lows_vec,:], dims=1)), vec(mean(firms_ie_col_SS[lows_vec,:], dims=1)), vec(mean(firms_y_col_SS[lows_vec,:], dims=1)), vec(mean(firms_Ei_col_SS[lows_vec,:], dims=1))
    firms_i_mean_low_IRF, firms_mp_mean_low_IRF, firms_ie_mean_low_IRF, firms_y_mean_low_IRF, firms_Ei_mean_low_IRF = vec(mean(firms_i_col_IRF[lows_vec,:], dims=1)), vec(mean(firms_mp_col_IRF[lows_vec,:], dims=1)), vec(mean(firms_ie_col_IRF[lows_vec,:], dims=1)), vec(mean(firms_y_col_IRF[lows_vec,:], dims=1)), vec(mean(firms_Ei_col_IRF[lows_vec,:], dims=1))
    # Also for log investment rates
    firms_li_mean_low_SS = vec(mean(log.(firms_i_col_SS[lows_vec,:]), dims=1))
    firms_li_mean_low_IRF = vec(mean(log.(firms_i_col_IRF[lows_vec,:]), dims=1))
    # Also for AC
    firms_liAC_mean_low_SS = vec(mean(log.(firms_iAC_col_SS[lows_vec,:]), dims=1))
    firms_liAC_mean_low_IRF = vec(mean(log.(firms_iAC_col_IRF[lows_vec,:]), dims=1))

    # Repeat for high subsample
    firms_i_mean_high_SS, firms_mp_mean_high_SS, firms_ie_mean_high_SS, firms_y_mean_high_SS, firms_Ei_mean_high_SS  = vec(mean(firms_i_col_SS[highs_vec,:], dims=1)), vec(mean(firms_mp_col_SS[highs_vec,:], dims=1)), vec(mean(firms_ie_col_SS[highs_vec,:], dims=1)), vec(mean(firms_y_col_SS[highs_vec,:], dims=1)), vec(mean(firms_Ei_col_SS[highs_vec,:], dims=1))
    firms_i_mean_high_IRF, firms_mp_mean_high_IRF, firms_ie_mean_high_IRF, firms_y_mean_high_IRF, firms_Ei_mean_high_IRF = vec(mean(firms_i_col_IRF[highs_vec,:], dims=1)), vec(mean(firms_mp_col_IRF[highs_vec,:], dims=1)), vec(mean(firms_ie_col_IRF[highs_vec,:], dims=1)), vec(mean(firms_y_col_IRF[highs_vec,:], dims=1)), vec(mean(firms_Ei_col_IRF[highs_vec,:], dims=1))
    # Also for log investment rates
    firms_li_mean_high_SS = vec(mean(log.(firms_i_col_SS[highs_vec,:]), dims=1))
    firms_li_mean_high_IRF = vec(mean(log.(firms_i_col_IRF[highs_vec,:]), dims=1))
    # Also for AC
    firms_liAC_mean_high_SS = vec(mean(log.(firms_iAC_col_SS[highs_vec,:]), dims=1))
    firms_liAC_mean_high_IRF = vec(mean(log.(firms_iAC_col_IRF[highs_vec,:]), dims=1))

    # Save all into collectors, given the current shock (SCALE RIGHT AWAY)
    firms_i_mean_IRFdev_col[is,:], firms_mp_mean_IRFdev_col[is,:], firms_ie_mean_IRFdev_col[is,:], firms_y_mean_IRFdev_col[is,:], firms_Ei_mean_IRFdev_col[is,:], firms_li_mean_IRFdev_col[is,:]  = (1/eps_scale)*(firms_i_mean_IRF - firms_i_mean_SS),  (1/eps_scale)*(firms_mp_mean_IRF - firms_mp_mean_SS),  (1/eps_scale)*(firms_ie_mean_IRF - firms_ie_mean_SS),  (1/eps_scale)*(firms_y_mean_IRF - firms_y_mean_SS),  (1/eps_scale)*(firms_Ei_mean_IRF - firms_Ei_mean_SS),  (1/eps_scale)*(firms_li_mean_IRF - firms_li_mean_SS)
    firms_i_mean_IRFdev_low_col[is,:], firms_mp_mean_IRFdev_low_col[is,:], firms_ie_mean_IRFdev_low_col[is,:], firms_y_mean_IRFdev_low_col[is,:], firms_Ei_mean_IRFdev_low_col[is,:], firms_li_mean_IRFdev_low_col[is,:]  = (1/eps_scale)*(firms_i_mean_low_IRF - firms_i_mean_low_SS),  (1/eps_scale)*(firms_mp_mean_low_IRF - firms_mp_mean_low_SS),  (1/eps_scale)*(firms_ie_mean_low_IRF - firms_ie_mean_low_SS),  (1/eps_scale)*(firms_y_mean_low_IRF - firms_y_mean_low_SS),  (1/eps_scale)*(firms_Ei_mean_low_IRF - firms_Ei_mean_low_SS),  (1/eps_scale)*(firms_li_mean_low_IRF - firms_li_mean_low_SS)
    firms_i_mean_IRFdev_high_col[is,:], firms_mp_mean_IRFdev_high_col[is,:], firms_ie_mean_IRFdev_high_col[is,:], firms_y_mean_IRFdev_high_col[is,:], firms_Ei_mean_IRFdev_high_col[is,:], firms_li_mean_IRFdev_high_col[is,:]  = (1/eps_scale)*(firms_i_mean_high_IRF - firms_i_mean_high_SS),  (1/eps_scale)*(firms_mp_mean_high_IRF - firms_mp_mean_high_SS),  (1/eps_scale)*(firms_ie_mean_high_IRF - firms_ie_mean_high_SS),  (1/eps_scale)*(firms_y_mean_high_IRF - firms_y_mean_high_SS),  (1/eps_scale)*(firms_Ei_mean_high_IRF - firms_Ei_mean_high_SS),  (1/eps_scale)*(firms_li_mean_high_IRF - firms_li_mean_high_SS)
    # Also for AC
    firms_liAC_mean_IRFdev_col[is,:] = (1/eps_scale)*(firms_liAC_mean_IRF - firms_liAC_mean_SS)
    firms_liAC_mean_IRFdev_low_col[is,:] = (1/eps_scale)*(firms_liAC_mean_low_IRF - firms_liAC_mean_low_SS)
    firms_liAC_mean_IRFdev_high_col[is,:] = (1/eps_scale)*(firms_liAC_mean_high_IRF - firms_liAC_mean_high_SS)
end

# Given the collections of log investment IRFs, compute the (weighted) average
w_matrix = repeat(sweights_vec,1,T)
firms_li_shmean_IRFdev      = sum(w_matrix .* firms_li_mean_IRFdev_col,dims=1)
firms_li_shmean_IRFdev_low  = sum(w_matrix .* firms_li_mean_IRFdev_low_col,dims=1)
firms_li_shmean_IRFdev_high = sum(w_matrix .* firms_li_mean_IRFdev_high_col,dims=1)
# Also for AC
firms_liAC_shmean_IRFdev      = sum(w_matrix .* firms_liAC_mean_IRFdev_col,dims=1)
firms_liAC_shmean_IRFdev_low  = sum(w_matrix .* firms_liAC_mean_IRFdev_low_col,dims=1)
firms_liAC_shmean_IRFdev_high = sum(w_matrix .* firms_liAC_mean_IRFdev_high_col,dims=1)

# Plot comparison with data
if !isdir("saved_figs/") mkdir("saved_figs/") end
if !isdir("saved_figs/_fig7") mkdir("saved_figs/_fig7") end

manl_ftsize = 13
ylim_man = (-2.9, 2.6)
plot_font = "helvetica"
default(fontfamily=plot_font, linewidth=2, framestyle=:box, label=nothing, grid=true)
# log investment rates
selr="high"
# Load empirical estimates
df_emp_betas= DataFrame(CSV.File("saved_figs/slevsbetas_vec_q_vals(linv_rat_realh)(ltobin_q)_IV_FE_INDT_FC_CLUftn_nwd.csv"))
df_emp_ci   = DataFrame(CSV.File("saved_figs/slevsbetas_vec_q_cis(linv_rat_realh)(ltobin_q)_IV_FE_INDT_FC_CLUftn_nwd.csv"))
emp_betas   = Array(df_emp_betas.V1)
emp_ci      = zeros(size(emp_betas,1),2)
emp_ci[:,1], emp_ci[:,2] = df_emp_ci.V1, df_emp_ci.V2
# Plot investment rate IRF-SS differences
Plots.plot(0:(T-1), [100*firms_liAC_shmean_IRFdev_high' emp_betas[1:T] emp_ci[1:T,1] emp_ci[1:T,2]], labels=["Model" "Data" "" ""], color=["green" "blue" "blue" "blue"], line = [:solid :solid :dash :dash], width=2.0, ylabel="Percent", xtickfontsize=manl_ftsize-1, ytickfontsize=manl_ftsize-1, yguidefontsize=manl_ftsize, xguidefontsize=manl_ftsize, legendfontsize=manl_ftsize, xticks=0:(T-1), xlabel="Quarters (h)", ylims=ylim_man, legend=:bottomleft)
Plots.hline!([0], color=:black, linestyle=:dash, labels="")
savefig("saved_figs/_fig7/ssdist-emp-"*selr*"_linvAC_sim_firm_shmeans_IRFdev_scaled_xirs.svg")

selr="low"
# Load empirical estimates
df_emp_betas= DataFrame(CSV.File("saved_figs/slevsbetas_ylevb_vec_q_vals(linv_rat_realh)(ltobin_q)_IV_FE_INDT_FC_CLUftn_nwd.csv"))
df_emp_ci   = DataFrame(CSV.File("saved_figs/slevsbetas_ylevb_vec_q_cis(linv_rat_realh)(ltobin_q)_IV_FE_INDT_FC_CLUftn_nwd.csv"))
emp_betas   = Array(df_emp_betas.V1)
emp_ci      = zeros(size(emp_betas,1),2)
emp_ci[:,1], emp_ci[:,2] = df_emp_ci.V1, df_emp_ci.V2
# Plot investment rate IRF-SS differences
Plots.plot(0:(T-1), [100*firms_liAC_shmean_IRFdev_low' emp_betas[1:T] emp_ci[1:T,1] emp_ci[1:T,2]], labels=["Model" "Data" "" ""], color=["green" "red" "red" "red"], line = [:solid :solid :dash :dash], width=2.0, ylabel="Percent", xtickfontsize=manl_ftsize-1, ytickfontsize=manl_ftsize-1, yguidefontsize=manl_ftsize, xguidefontsize=manl_ftsize, legendfontsize=manl_ftsize, xticks=0:(T-1), xlabel="Quarters (h)", ylims=ylim_man, legend=:bottomleft)
Plots.hline!([0], color=:black, linestyle=:dash, labels="")
savefig("saved_figs/_fig7/ssdist-emp-"*selr*"_linvAC_sim_firm_shmeans_IRFdev_scaled_xirs.svg")


###
# If requested, also compute calibration moments (As seen in Table 1)
###
if false
    using LinearAlgebra, SparseArrays

    # Compute SS policies on fine grid
    vf_sol1, = solve_valfunc_v(c_sol, modl.psi, modl.rm, sspace.sf0, modl, sspace, opts, false)
    vf_sol0, = solve_valfunc_v0(c_sol, modl.psi, modl.rm, sspace.sf0, modl, sspace, opts, false)
    # Compute cost cutoffs
    xi_vecf  = vf_sol1.val - vf_sol0.val
    # Probabilities of issuing
    Ei_prob  = xi_cdfun(xi_vecf, modl, sspace, opts)

    # Create function for computing quantiles of y, given discrete grid and distribution
    function discr_percs(ygrid, L, percs)
        isort_ygrid = sortperm(ygrid)
        ygrid_sort  = ygrid[isort_ygrid]
        L_sort      = L[isort_ygrid]
        # Create cumulative distribution
        L_ysort_cum = cumsum(L_sort)
        y_perc_vals = zeros(length(percs))
        for ip=1:length(percs)
            # Determine a specific percentile
            y_perc         = percs[ip]
            # The percentile value
            y_perc_i       = findfirst(L_ysort_cum .> y_perc)
            if y_perc_i == 1
                y_perc_val = ygrid_sort[1]
            else
                y_perc_val     = ygrid_sort[y_perc_i-1] + (ygrid_sort[y_perc_i]-ygrid_sort[y_perc_i-1]) * (y_perc-L_ysort_cum[y_perc_i-1])/(L_ysort_cum[y_perc_i]-L_ysort_cum[y_perc_i-1])
            end
            y_perc_vals[ip]= y_perc_val
        end
        return y_perc_vals
    end

    # Construct marginal distribution of cash to assets
    Nd  = size(L_sol,1)
    Lm  = kron(ones(1,sspace.nf[3]), I(sspace.nf[2])) * L_sol
    # 1. Median M/A
    med_m   = discr_percs(sspace.mgridf, Lm, 0.5)[1] # Median M/K
    med_liq = med_m/(1.0+med_m)                      # Median M/A
    # 2. and 3. Median I/K for low and high M/A
    # Given the median M/K, compute high/low conditional distributions
    L_lowm_long = vcat(Ei_prob .* L_sol .* (sspace.sf0[:,1] .<= med_m), (ones(Nd) - Ei_prob).* L_sol .* (sspace.sf0[:,1] .<= med_m))
    L_highm_long = vcat(Ei_prob .* L_sol .* (sspace.sf0[:,1] .> med_m), (ones(Nd) - Ei_prob).* L_sol .* (sspace.sf0[:,1] .> med_m))
    L_lowm_long = L_lowm_long/sum(L_lowm_long)
    L_highm_long= L_highm_long/sum(L_highm_long)
    # Compute medians
    i_pol_long  = vcat(vf_sol1.i, vf_sol0.i)
    med_i_lowm   = discr_percs(i_pol_long, L_lowm_long, 0.5)[1]
    med_i_highm  = discr_percs(i_pol_long, L_highm_long, 0.5)[1]
    # Including adjustment cost
    med_iAC_lowm   = med_i_lowm + modl.nu*((med_i_lowm - modl.delta_0)^2)
    med_iAC_highm  = med_i_highm + modl.nu*((med_i_highm - modl.delta_0)^2)
    # 4. Frequency of issuance
    freq_eqiss  = L_sol' * (Ei_prob.*(vf_sol1.ie .> 0.05/4)) + L_sol' * ((ones(Nd)-Ei_prob).*(vf_sol0.ie .> 0.05/4))
    # 5. Mean equity issuance, conditional on issuing
    L_Epos1    = (Ei_prob.* L_sol) .* (vf_sol1.ie .> 0.05/4)
    L_Epos0    = ((ones(Nd)-Ei_prob) .* L_sol) .* (vf_sol0.ie .> 0.05/4)
    avg_EdA_Epos= (L_Epos1' * (vf_sol1.ie*modl.psi ./ (ones(size(sspace.sf0[:,1],1)) + sspace.sf0[:,1])) + L_Epos0' * (vf_sol0.ie*modl.psi ./ (ones(size(sspace.sf0[:,1],1)) + sspace.sf0[:,1]))) / sum(L_Epos1 + L_Epos0)

    println("### Calibration moments ###")
    println("median(M/A): $(round(med_liq,digits=3))")
    println("median(I/K)| low M/A: $(round(med_iAC_lowm, digits=3))")
    println("median(I/K)| high M/A: $(round(med_iAC_highm, digits=3))")
    println("frequency of issuance: $(round(freq_eqiss,digits=3))")
    println("mean(e)| issuance: $(round(avg_EdA_Epos,digits=3))")
end
####
